knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE #

if (get_bayesfactor) {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = paste0('_AsPreregisteredAPIM_', as.character(iterations))
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1649 0.0000 5.0000 275

Prepare for APIM

df_double$is_A <- ifelse(df_double$pa_obj_self_cb >= df_double$pa_obj_partner_cb, 0, 1)
df_double$is_B <- ifelse(df_double$pa_obj_self_cb < df_double$pa_obj_partner_cb, 0, 1)

Modelling

Self-Reported MVPA

Persuasion

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID),
  
  hu = ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
  , brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub_persuasion <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_persuasion", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
## Start sampling
if (check_models) {
  check_brms(pa_sub_persuasion, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub_persuasion, integer = TRUE, outliers_type = 'bootstrap')
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                        Term  VIF   VIF 95% CI Increased SE Tolerance
##                        is_A 2.24 [2.14, 2.33]         1.49      0.45
##                        is_B 1.06 [1.03, 1.10]         1.03      0.95
##     is_A:persuasion_self_cw 1.15 [1.12, 1.19]         1.07      0.87
##     is_B:persuasion_self_cw 1.14 [1.11, 1.18]         1.07      0.88
##  is_A:persuasion_partner_cw 1.04 [1.02, 1.09]         1.02      0.96
##  is_B:persuasion_partner_cw 2.71 [2.60, 2.84]         1.65      0.37
##     is_A:persuasion_self_cb 2.86 [2.74, 2.99]         1.69      0.35
##  Tolerance 95% CI
##      [0.43, 0.47]
##      [0.91, 0.97]
##      [0.84, 0.89]
##      [0.85, 0.90]
##      [0.92, 0.98]
##      [0.35, 0.38]
##      [0.33, 0.37]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         94%
##     lognormal          6%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##  neg. binomial (zero-infl.)         78%
##               beta-binomial         16%
##                   lognormal          3%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 7, observations = 3742, p-value =
## 0.84
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0005344735 0.0029396045
## sample estimates:
## outlier frequency (expected: 0.001574024585783 ) 
##                                      0.001870657
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_persuasion$data$pa_sub[pa_sub_persuasion$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub_persuasion, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub_persuasion, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
is_A 0.74 0.14 [ 0.51, 1.07] 0.944 [0.84, 1.20] 0.257 0.260 Moderate Evidence for Null 1.000 17468 24772 46.58*** 3.50 [40.07, 53.85] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 13642 22685
is_A:day 0.97 0.18 [ 0.67, 1.40] 0.570 [0.84, 1.20] 0.660 0.095 Strong Evidence for Null 1.000 44612 29960 1.09 0.09 [ 0.92, 1.29] 0.845 [0.92, 1.08] 0.436 0.012 Very Strong Evidence for Null 1.000 46827 29347
is_A:persuasion_partner_cb 0.75 0.30 [ 0.34, 1.64] 0.767 [0.84, 1.20] 0.271 0.254 Moderate Evidence for Null 1.000 14088 21142 1.02 0.17 [ 0.74, 1.42] 0.561 [0.92, 1.08] 0.365 0.018 Very Strong Evidence for Null 1.000 11377 19396
is_A:persuasion_partner_cw 1.58*** 0.16 [ 1.30, 1.98] 1.000 [0.84, 1.20] 0.004 >100 Overwhelming Evidence 1.000 23821 26789 1.02 0.03 [ 0.96, 1.08] 0.712 [0.92, 1.08] 0.980 0.008 Very Strong Evidence for Null 1.000 20463 25811
is_A:persuasion_self_cb 0.81 0.36 [ 0.33, 1.97] 0.683 [0.84, 1.20] 0.281 0.251 Moderate Evidence for Null 1.000 13340 20895 1.21 0.21 [ 0.84, 1.72] 0.853 [0.92, 1.08] 0.199 0.029 Very Strong Evidence for Null 1.001 10250 19030
is_A:persuasion_self_cw 1.67*** 0.16 [ 1.41, 2.10] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 22586 23823 1.01 0.04 [ 0.94, 1.10] 0.643 [0.92, 1.08] 0.936 0.011 Very Strong Evidence for Null 1.000 15189 22445
is_B 0.92 0.18 [ 0.63, 1.37] 0.656 [0.84, 1.20] 0.597 0.084 Strong Evidence for Null 1.000 16516 22811 49.76*** 3.97 [42.49, 58.38] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 11733 19760
is_B:day 0.87 0.16 [ 0.61, 1.25] 0.777 [0.84, 1.20] 0.543 0.123 Moderate Evidence for Null 1.000 41627 29767 0.88 0.07 [ 0.75, 1.03] 0.943 [0.92, 1.08] 0.257 0.024 Very Strong Evidence for Null 1.000 50533 30069
is_B:persuasion_partner_cb 1.01 0.48 [ 0.39, 2.62] 0.510 [0.84, 1.20] 0.293 0.235 Moderate Evidence for Null 1.000 12851 19993 1.03 0.19 [ 0.72, 1.48] 0.568 [0.92, 1.08] 0.330 0.017 Very Strong Evidence for Null 1.000 11737 19703
is_B:persuasion_partner_cw 1.49*** 0.14 [ 1.24, 1.85] 1.000 [0.84, 1.20] 0.010 86.845 Very Strong Evidence 1.000 26845 24277 1.03 0.03 [ 0.97, 1.10] 0.858 [0.92, 1.08] 0.935 0.014 Very Strong Evidence for Null 1.000 20439 27855
is_B:persuasion_self_cb 0.70 0.29 [ 0.30, 1.62] 0.805 [0.84, 1.20] 0.233 0.304 Weak Evidence for Null 1.000 13544 21173 1.13 0.18 [ 0.82, 1.57] 0.767 [0.92, 1.08] 0.289 0.023 Very Strong Evidence for Null 1.000 12920 20807
is_B:persuasion_self_cw 1.68*** 0.17 [ 1.39, 2.10] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 22677 24944 1.04 0.03 [ 0.98, 1.10] 0.889 [0.92, 1.08] 0.930 0.016 Very Strong Evidence for Null 1.000 17302 24439
sd(is_A) 0.94 0.12 [0.73, 1.23] NA NA NA NA NA 1.000 20183 25071 0.33 0.05 [0.25, 0.44] NA NA NA NA NA 1.000 13315 23621
sd(is_A:persuasion_partner_cw) 0.41 0.12 [0.21, 0.69] NA NA NA NA NA 1.000 18656 21846 0.06 0.03 [0.01, 0.13] NA NA NA NA NA 1.000 12477 10932
sd(is_A:persuasion_self_cw) 0.29 0.13 [0.05, 0.58] NA NA NA NA NA 1.000 9601 8514 0.17 0.03 [0.11, 0.24] NA NA NA NA NA 1.000 20343 28866
sd(is_B) 1.01 0.13 [0.79, 1.32] NA NA NA NA NA 1.000 18297 24221 0.39 0.05 [0.30, 0.50] NA NA NA NA NA 1.000 13771 21334
sd(is_B:persuasion_partner_cw) 0.33 0.14 [0.06, 0.63] NA NA NA NA NA 1.000 10537 8117 0.10 0.03 [0.05, 0.16] NA NA NA NA NA 1.000 20368 18826
sd(is_B:persuasion_self_cw) 0.38 0.10 [0.19, 0.62] NA NA NA NA NA 1.000 21344 22705 0.09 0.03 [0.05, 0.15] NA NA NA NA NA 1.000 18840 17891
sigma NA NA NA NA NA NA NA NA NA NA NA 0.67 0.01 [0.65, 0.70] NA NA NA NA NA 1.000 47834 30606

Hurdle part of the model on the left, non-zero part towards the right side of the table

Pressure

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID),
  
  hu = ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
  , brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub_pressure <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_pressure", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
## Start sampling
if (check_models) {
  check_brms(pa_sub_pressure, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub_pressure, integer = TRUE, outliers_type = 'bootstrap')
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                      Term  VIF   VIF 95% CI Increased SE Tolerance
##                      is_A 2.23 [2.14, 2.33]         1.49      0.45
##                      is_B 1.05 [1.02, 1.09]         1.02      0.96
##     is_A:pressure_self_cw 1.01 [1.00, 1.14]         1.01      0.99
##     is_B:pressure_self_cw 1.01 [1.00, 1.12]         1.01      0.99
##  is_A:pressure_partner_cw 1.01 [1.00, 1.15]         1.01      0.99
##  is_B:pressure_partner_cw 3.49 [3.33, 3.66]         1.87      0.29
##     is_A:pressure_self_cb 4.30 [4.10, 4.52]         2.07      0.23
##  Tolerance 95% CI
##      [0.43, 0.47]
##      [0.92, 0.98]
##      [0.88, 1.00]
##      [0.89, 1.00]
##      [0.87, 1.00]
##      [0.27, 0.30]
##      [0.22, 0.24]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         94%
##     lognormal          6%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##  neg. binomial (zero-infl.)         78%
##               beta-binomial         16%
##                   lognormal          3%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 11 observations with a pareto_k > 0.7 in model 'model'.
## We recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 6, observations = 3736, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0005353319 0.0030848501
## sample estimates:
## outlier frequency (expected: 0.00157922912205567 ) 
##                                        0.001605996
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_pressure$data$pa_sub[pa_sub_pressure$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub_pressure, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub_pressure, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
is_A 0.87 0.15 [ 0.61, 1.23] 0.783 [0.84, 1.20] 0.555 0.096 Strong Evidence for Null 1.000 14549 20069 48.25*** 3.51 [41.79, 55.67] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 10797 19767
is_A:day 0.67* 0.12 [ 0.48, 0.95] 0.989 [0.84, 1.20] 0.107 1.172 Weak Evidence 1.000 37555 29688 1.09 0.09 [ 0.92, 1.28] 0.835 [0.92, 1.08] 0.452 0.011 Very Strong Evidence for Null 1.000 48844 30086
is_A:pressure_partner_cb 0.63 0.39 [ 0.18, 2.10] 0.772 [0.84, 1.20] 0.174 0.412 Weak Evidence for Null 1.000 14001 22562 0.94 0.26 [ 0.54, 1.64] 0.582 [0.92, 1.08] 0.218 0.016 Very Strong Evidence for Null 1.000 12801 20013
is_A:pressure_partner_cw 2.44*** 0.57 [ 1.57, 4.32] 1.000 [0.84, 1.20] 0.001 >100 Overwhelming Evidence 1.000 27317 18821 0.93 0.05 [ 0.84, 1.05] 0.880 [0.92, 1.08] 0.567 0.012 Very Strong Evidence for Null 1.000 29892 25384
is_A:pressure_self_cb 0.57 0.44 [ 0.13, 2.62] 0.766 [0.84, 1.20] 0.143 0.502 Weak Evidence for Null 1.000 14022 21390 1.56 0.53 [ 0.80, 3.08] 0.908 [0.92, 1.08] 0.075 0.038 Strong Evidence for Null 1.000 14697 21985
is_A:pressure_self_cw 1.86* 0.62 [ 1.02, 4.68] 0.978 [0.84, 1.20] 0.060 1.518 Weak Evidence 1.000 19579 18590 0.85* 0.06 [ 0.72, 1.00] 0.978 [0.92, 1.08] 0.136 0.082 Strong Evidence for Null 1.000 23058 21388
is_B 1.14 0.21 [ 0.80, 1.63] 0.763 [0.84, 1.20] 0.560 0.090 Strong Evidence for Null 1.000 13997 21736 52.20*** 4.30 [44.46, 61.41] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.001 9866 16264
is_B:day 0.56*** 0.10 [ 0.40, 0.79] 1.000 [0.84, 1.20] 0.011 20.363 Strong Evidence 1.000 36258 29993 0.83* 0.07 [ 0.71, 0.97] 0.991 [0.92, 1.08] 0.085 0.109 Moderate Evidence for Null 1.000 48596 30206
is_B:pressure_partner_cb 1.13 0.88 [ 0.24, 5.28] 0.562 [0.84, 1.20] 0.183 0.393 Weak Evidence for Null 1.000 12957 20347 1.10 0.43 [ 0.49, 2.47] 0.593 [0.92, 1.08] 0.155 0.018 Very Strong Evidence for Null 1.001 12994 19397
is_B:pressure_partner_cw 3.28* 2.01 [ 1.24, 16.35] 0.990 [0.84, 1.20] 0.017 4.367 Moderate Evidence 1.000 17196 17609 0.99 0.07 [ 0.82, 1.13] 0.575 [0.92, 1.08] 0.709 0.008 Very Strong Evidence for Null 1.000 21016 19395
is_B:pressure_self_cb 0.48 0.30 [ 0.14, 1.67] 0.882 [0.84, 1.20] 0.111 0.645 Weak Evidence for Null 1.000 12815 19815 1.08 0.34 [ 0.57, 2.05] 0.593 [0.92, 1.08] 0.194 0.019 Very Strong Evidence for Null 1.001 12751 19790
is_B:pressure_self_cw 1.29 0.34 [ 0.71, 2.66] 0.832 [0.84, 1.20] 0.326 0.227 Moderate Evidence for Null 1.000 20759 17702 0.93 0.06 [ 0.82, 1.08] 0.842 [0.92, 1.08] 0.545 0.012 Very Strong Evidence for Null 1.000 23751 24390
sd(is_A) 0.88 0.12 [0.68, 1.15] NA NA NA NA NA 1.000 17349 24619 0.32 0.05 [0.24, 0.44] NA NA NA NA NA 1.000 12153 20889
sd(is_A:pressure_partner_cw) 0.31 0.28 [0.01, 1.20] NA NA NA NA NA 1.000 15331 16422 0.06 0.05 [0.00, 0.21] NA NA NA NA NA 1.000 17943 17235
sd(is_A:pressure_self_cw) 0.72 0.53 [0.06, 2.40] NA NA NA NA NA 1.000 11421 12927 0.10 0.09 [0.00, 0.38] NA NA NA NA NA 1.000 12722 15622
sd(is_B) 0.91 0.12 [0.70, 1.19] NA NA NA NA NA 1.000 18280 24673 0.40 0.05 [0.31, 0.53] NA NA NA NA NA 1.000 11041 18555
sd(is_B:pressure_partner_cw) 1.65 0.85 [0.31, 3.83] NA NA NA NA NA 1.000 12805 11132 0.09 0.08 [0.00, 0.37] NA NA NA NA NA 1.000 13239 14438
sd(is_B:pressure_self_cw) 0.74 0.49 [0.06, 2.12] NA NA NA NA NA 1.000 10281 10816 0.09 0.07 [0.00, 0.28] NA NA NA NA NA 1.000 14609 13572
sigma NA NA NA NA NA NA NA NA NA NA NA 0.69 0.01 [0.67, 0.72] NA NA NA NA NA 1.000 43470 29205

Hurdle part of the model on the left, non-zero part towards the right side of the table

Pushing

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    barriers_self_cw + barriers_self_cb + 
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID),
  
  hu = ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
  , brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub_pushing <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_pushing", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
## Start sampling
if (check_models) {
  check_brms(pa_sub_pushing, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub_pushing, integer = TRUE, outliers_type = 'bootstrap')
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                     Term  VIF   VIF 95% CI Increased SE Tolerance
##                     is_A 1.75 [1.68, 1.84]         1.32      0.57
##                     is_B 1.08 [1.05, 1.12]         1.04      0.93
##         barriers_self_cw 1.04 [1.02, 1.10]         1.02      0.96
##         barriers_self_cb 1.04 [1.01, 1.09]         1.02      0.97
##     is_A:pushing_self_cw 1.07 [1.04, 1.12]         1.04      0.93
##     is_B:pushing_self_cw 1.22 [1.18, 1.28]         1.11      0.82
##  is_A:pushing_partner_cw 1.24 [1.20, 1.30]         1.12      0.80
##  is_B:pushing_partner_cw 1.19 [1.15, 1.24]         1.09      0.84
##     is_A:pushing_self_cb 1.27 [1.23, 1.33]         1.13      0.79
##  Tolerance 95% CI
##      [0.54, 0.60]
##      [0.89, 0.95]
##      [0.91, 0.98]
##      [0.92, 0.99]
##      [0.89, 0.96]
##      [0.78, 0.85]
##      [0.77, 0.83]
##      [0.81, 0.87]
##      [0.75, 0.82]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         91%
##     lognormal          9%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##  neg. binomial (zero-infl.)         84%
##               beta-binomial          9%
##                   lognormal          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 12 observations with a pareto_k > 0.7 in model 'model'.
## We recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 5, observations = 1776, p-value =
## 0.36
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.00423705
## sample estimates:
## outlier frequency (expected: 0.00161036036036036 ) 
##                                        0.002815315
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_pushing$data$pa_sub[pa_sub_pushing$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub_pushing, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub_pushing, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
barriers_self_cb NA NA NA NA NA NA NA NA NA NA NA 0.89 0.13 [ 0.68, 1.18] 0.791 [0.93, 1.08] 0.305 0.026 Very Strong Evidence for Null 1.000 15389 22099
barriers_self_cw NA NA NA NA NA NA NA NA NA NA NA 1.06 0.04 [ 0.99, 1.13] 0.943 [0.93, 1.08] 0.739 0.032 Strong Evidence for Null 1.000 68293 29503
is_A 2.63*** 0.70 [ 1.56, 4.51] 1.000 [0.84, 1.20] 0.002 48.662 Very Strong Evidence 1.000 29400 28949 53.39*** 4.54 [45.24, 63.19] 1.000 [0.93, 1.08] 0.000 >100 Overwhelming Evidence 1.000 18788 26280
is_A:day 0.82 0.24 [ 0.46, 1.45] 0.751 [0.84, 1.20] 0.369 0.187 Moderate Evidence for Null 1.000 58113 28964 1.08 0.10 [ 0.90, 1.30] 0.793 [0.93, 1.08] 0.441 0.011 Very Strong Evidence for Null 1.000 66370 30302
is_A:pushing_partner_cb 0.91 0.88 [ 0.13, 6.16] 0.540 [0.84, 1.20] 0.147 0.491 Weak Evidence for Null 1.000 30887 30885 1.10 0.40 [ 0.52, 2.26] 0.597 [0.93, 1.08] 0.159 0.020 Very Strong Evidence for Null 1.000 15843 23109
is_A:pushing_partner_cw 1.95** 0.45 [ 1.30, 3.48] 0.999 [0.84, 1.20] 0.009 14.509 Strong Evidence 1.000 32690 26682 0.96 0.03 [ 0.89, 1.03] 0.886 [0.93, 1.08] 0.837 0.017 Very Strong Evidence for Null 1.000 47103 27861
is_A:pushing_self_cb 0.39 0.35 [ 0.07, 2.31] 0.852 [0.84, 1.20] 0.093 0.769 Weak Evidence for Null 1.000 28601 30135 1.29 0.42 [ 0.67, 2.49] 0.785 [0.93, 1.08] 0.136 0.025 Very Strong Evidence for Null 1.000 18803 26620
is_A:pushing_self_cw 1.91 0.76 [ 0.96, 5.60] 0.967 [0.84, 1.20] 0.081 0.998 Weak Evidence for Null 1.000 25693 21815 1.00 0.06 [ 0.89, 1.12] 0.523 [0.93, 1.08] 0.814 0.012 Very Strong Evidence for Null 1.000 24382 27449
is_B 4.02*** 1.18 [ 2.24, 7.28] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 26916 28674 52.54*** 4.45 [44.44, 62.18] 1.000 [0.93, 1.08] 0.000 >100 Overwhelming Evidence 1.000 20116 25664
is_B:day 0.65 0.19 [ 0.36, 1.15] 0.929 [0.84, 1.20] 0.177 0.427 Weak Evidence for Null 1.000 63463 30099 0.94 0.09 [ 0.79, 1.12] 0.746 [0.93, 1.08] 0.507 0.010 Very Strong Evidence for Null 1.000 74905 30566
is_B:pushing_partner_cb 0.81 0.81 [ 0.11, 5.67] 0.582 [0.84, 1.20] 0.140 0.508 Weak Evidence for Null 1.000 29944 29697 1.04 0.34 [ 0.53, 2.01] 0.543 [0.93, 1.08] 0.187 0.020 Very Strong Evidence for Null 1.000 19801 25578
is_B:pushing_partner_cw 2.28* 1.11 [ 1.04, 8.55] 0.979 [0.84, 1.20] 0.044 1.868 Weak Evidence 1.000 21839 21728 1.02 0.06 [ 0.91, 1.15] 0.658 [0.93, 1.08] 0.797 0.012 Very Strong Evidence for Null 1.000 27666 30385
is_B:pushing_self_cb 0.41 0.42 [ 0.05, 3.18] 0.810 [0.84, 1.20] 0.096 0.757 Weak Evidence for Null 1.000 31398 29450 1.50 0.55 [ 0.71, 3.13] 0.859 [0.93, 1.08] 0.092 0.038 Strong Evidence for Null 1.000 17644 25436
is_B:pushing_self_cw 1.16 0.25 [ 0.73, 1.91] 0.753 [0.84, 1.20] 0.491 0.137 Moderate Evidence for Null 1.000 29855 26726 0.95 0.03 [ 0.88, 1.01] 0.939 [0.93, 1.08] 0.726 0.030 Very Strong Evidence for Null 1.000 43301 31543
sd(is_A) 1.18 0.20 [0.85, 1.64] NA NA NA NA NA 1.000 27290 31419 0.36 0.06 [0.26, 0.50] NA NA NA NA NA 1.000 16350 24786
sd(is_A:pushing_partner_cw) 0.65 0.33 [0.09, 1.45] NA NA NA NA NA 1.000 14902 13091 0.04 0.04 [0.00, 0.14] NA NA NA NA NA 1.000 20863 21692
sd(is_A:pushing_self_cw) 1.51 0.59 [0.49, 3.01] NA NA NA NA NA 1.000 13003 10459 0.18 0.06 [0.05, 0.32] NA NA NA NA NA 1.000 12931 9839
sd(is_B) 1.37 0.23 [0.99, 1.92] NA NA NA NA NA 1.000 25542 30110 0.38 0.06 [0.28, 0.52] NA NA NA NA NA 1.000 16644 25124
sd(is_B:pushing_partner_cw) 1.72 0.82 [0.24, 3.65] NA NA NA NA NA 1.000 10498 8114 0.17 0.07 [0.03, 0.32] NA NA NA NA NA 1.000 12628 10344
sd(is_B:pushing_self_cw) 0.85 0.33 [0.29, 1.72] NA NA NA NA NA 1.000 14504 13114 0.05 0.04 [0.00, 0.15] NA NA NA NA NA 1.000 19509 21091
sigma NA NA NA NA NA NA NA NA NA NA NA 0.67 0.01 [0.64, 0.70] NA NA NA NA NA 1.000 57411 29137

Hurdle part of the model on the left, non-zero part towards the right side of the table

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Persuasion

formula <- bf(
  pa_obj ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log_persuasion <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_persuasion", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
  check_brms(pa_obj_log_persuasion, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log_persuasion, integer = TRUE, outliers_type = 'bootstrap')
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                        Term  VIF   VIF 95% CI Increased SE Tolerance
##                        is_A 1.84 [1.76, 1.93]         1.36      0.54
##                        is_B 1.08 [1.05, 1.12]         1.04      0.93
##     is_A:persuasion_self_cw 1.13 [1.10, 1.18]         1.06      0.88
##     is_B:persuasion_self_cw 1.09 [1.06, 1.13]         1.04      0.92
##  is_A:persuasion_partner_cw 1.04 [1.02, 1.09]         1.02      0.96
##  is_B:persuasion_partner_cw 3.53 [3.35, 3.73]         1.88      0.28
##     is_A:persuasion_self_cb 2.94 [2.79, 3.10]         1.71      0.34
##     is_B:persuasion_self_cb 3.69 [3.50, 3.90]         1.92      0.27
##  is_A:persuasion_partner_cb 2.95 [2.80, 3.11]         1.72      0.34
##  Tolerance 95% CI
##      [0.52, 0.57]
##      [0.89, 0.95]
##      [0.85, 0.91]
##      [0.88, 0.95]
##      [0.91, 0.98]
##      [0.27, 0.30]
##      [0.32, 0.36]
##      [0.26, 0.29]
##      [0.32, 0.36]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         88%
##     lognormal          9%
##       weibull          3%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##                   lognormal         72%
##                     tweedie         16%
##  neg. binomial (zero-infl.)          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 33, observations = 3343, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.001338618 0.005085253
## sample estimates:
## outlier frequency (expected: 0.00303021238408615 ) 
##                                        0.009871373
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_persuasion$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log_persuasion, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log_persuasion, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
summary_pa_obj %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
is_A 101.53*** 6.28 [ 89.80, 114.88] 1.000 [0.94, 1.07] 0.000 >100 Overwhelming Evidence 1.000 9441 16504
is_B 135.95*** 7.32 [122.13, 151.50] 1.000 [0.94, 1.07] 0.000 >100 Overwhelming Evidence 1.000 13358 20795
is_A:persuasion_self_cw 1.05* 0.02 [ 1.01, 1.09] 0.989 [0.94, 1.07] 0.839 0.083 Strong Evidence for Null 1.000 36995 32045
is_B:persuasion_self_cw 1.01 0.02 [ 0.97, 1.05] 0.659 [0.94, 1.07] 0.998 0.005 Very Strong Evidence for Null 1.000 35355 33024
is_A:persuasion_partner_cw 1.02 0.02 [ 0.98, 1.06] 0.792 [0.94, 1.07] 0.988 0.007 Very Strong Evidence for Null 1.000 36380 30708
is_B:persuasion_partner_cw 1.03 0.02 [ 0.99, 1.07] 0.909 [0.94, 1.07] 0.969 0.014 Very Strong Evidence for Null 1.000 40799 30254
is_A:persuasion_self_cb 1.25 0.21 [ 0.89, 1.76] 0.901 [0.94, 1.07] 0.128 0.033 Strong Evidence for Null 1.000 10402 17765
is_B:persuasion_self_cb 0.89 0.11 [ 0.70, 1.14] 0.824 [0.94, 1.07] 0.263 0.019 Very Strong Evidence for Null 1.000 14430 21279
is_A:persuasion_partner_cb 0.89 0.13 [ 0.66, 1.20] 0.787 [0.94, 1.07] 0.249 0.019 Very Strong Evidence for Null 1.000 11805 19041
is_B:persuasion_partner_cb 1.29 0.19 [ 0.96, 1.73] 0.957 [0.94, 1.07] 0.081 0.062 Strong Evidence for Null 1.000 13364 21848
is_A:day 0.98 0.05 [ 0.90, 1.08] 0.628 [0.94, 1.07] 0.798 0.004 Very Strong Evidence for Null 1.000 87762 28942
is_B:day 0.95 0.04 [ 0.87, 1.04] 0.859 [0.94, 1.07] 0.612 0.007 Very Strong Evidence for Null 1.000 93845 28595
is_A:weartime_self_cw 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.07] 1.000 11.432 Strong Evidence 1.000 79311 30606
is_B:weartime_self_cw 1.00* 0.00 [ 1.00, 1.00] 0.983 [0.94, 1.07] 1.000 0.039 Strong Evidence for Null 1.000 80060 29963
is_A:weartime_self_cb 1.00 0.00 [ 1.00, 1.00] 0.712 [0.94, 1.07] 1.000 0.014 Very Strong Evidence for Null 1.000 12377 18753
is_B:weartime_self_cb 1.00 0.00 [ 1.00, 1.00] 0.584 [0.94, 1.07] 1.000 0.011 Very Strong Evidence for Null 1.000 15565 23274
sd(is_A) 0.34 0.04 [0.26, 0.44] NA NA NA NA NA 1.000 12417 22026
sd(is_B) 0.28 0.04 [0.22, 0.37] NA NA NA NA NA 1.000 14817 22063
sd(is_A:persuasion_self_cw) 0.06 0.02 [0.02, 0.11] NA NA NA NA NA 1.000 14879 11539
sd(is_B:persuasion_self_cw) 0.06 0.02 [0.02, 0.11] NA NA NA NA NA 1.000 17283 14227
sd(is_A:persuasion_partner_cw) 0.07 0.02 [0.03, 0.12] NA NA NA NA NA 1.000 15202 15518
sd(is_B:persuasion_partner_cw) 0.07 0.03 [0.01, 0.13] NA NA NA NA NA 1.000 11803 13070
sigma 0.55 0.01 [0.53, 0.56] NA NA NA NA NA 1.000 73910 28826

Pressure

formula <- bf(
  pa_obj ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log_pressure <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_pressure", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
  check_brms(pa_obj_log_pressure, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log_pressure, integer = TRUE, outliers_type = 'bootstrap')
}
## Model has no intercept. VIFs may not be sensible.
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model
##   without interaction terms.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                      Term  VIF     VIF 95% CI Increased SE Tolerance
##                      is_A 1.86 [ 1.78,  1.95]         1.36      0.54
##                      is_B 1.02 [ 1.01,  1.09]         1.01      0.98
##     is_A:pressure_self_cw 1.01 [ 1.00,  1.13]         1.01      0.99
##     is_B:pressure_self_cw 1.02 [ 1.00,  1.11]         1.01      0.98
##  is_A:pressure_partner_cw 1.04 [ 1.01,  1.09]         1.02      0.96
##  Tolerance 95% CI
##      [0.51, 0.56]
##      [0.91, 0.99]
##      [0.88, 1.00]
##      [0.90, 1.00]
##      [0.92, 0.99]
## 
## Moderate Correlation
## 
##                      Term  VIF     VIF 95% CI Increased SE Tolerance
##     is_A:pressure_self_cb 8.49 [ 8.01,  9.01]         2.91      0.12
##  is_A:pressure_partner_cb 8.81 [ 8.30,  9.34]         2.97      0.11
##  Tolerance 95% CI
##      [0.11, 0.12]
##      [0.11, 0.12]
## 
## High Correlation
## 
##                      Term   VIF     VIF 95% CI Increased SE Tolerance
##  is_B:pressure_partner_cw 10.71 [10.10, 11.37]         3.27      0.09
##     is_B:pressure_self_cb 10.70 [10.08, 11.35]         3.27      0.09
##  Tolerance 95% CI
##      [0.09, 0.10]
##      [0.09, 0.10]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         88%
##     lognormal          9%
##       weibull          3%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##                   lognormal         72%
##                     tweedie         16%
##  neg. binomial (zero-infl.)          9%
## 
## Divergences:
## 2 of 40000 iterations ended with a divergence (0.005%).
## Try increasing 'adapt_delta' to remove the divergences.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 34, observations = 3337, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.001640695 0.004794726
## sample estimates:
## outlier frequency (expected: 0.00308660473479173 ) 
##                                         0.01018879
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_pressure$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log_pressure, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Warning: There were 2 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Model has no intercept. VIFs may not be sensible.
## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model
##   without interaction terms.
## Warning in compute_rope(range = rope_range): Collinearity detected. Some
## VIFs are > 10. This may invalidate ROPE inferences!

Check for Multicollinearity

Low Correlation

                 Term  VIF     VIF 95% CI Increased SE Tolerance
                 is_A 1.86 [ 1.78,  1.95]         1.36      0.54
                 is_B 1.02 [ 1.01,  1.09]         1.01      0.98
is_A:pressure_self_cw 1.01 [ 1.00,  1.13]         1.01      0.99
is_B:pressure_self_cw 1.02 [ 1.00,  1.11]         1.01      0.98

is_A:pressure_partner_cw 1.04 [ 1.01, 1.09] 1.02 0.96 Tolerance 95% CI [0.51, 0.56] [0.91, 0.99] [0.88, 1.00] [0.90, 1.00] [0.92, 0.99]

Moderate Correlation

                 Term  VIF     VIF 95% CI Increased SE Tolerance
is_A:pressure_self_cb 8.49 [ 8.01,  9.01]         2.91      0.12

is_A:pressure_partner_cb 8.81 [ 8.30, 9.34] 2.97 0.11 Tolerance 95% CI [0.11, 0.12] [0.11, 0.12]

High Correlation

                 Term   VIF     VIF 95% CI Increased SE Tolerance

is_B:pressure_partner_cw 10.71 [10.10, 11.37] 3.27 0.09 is_B:pressure_self_cb 10.70 [10.08, 11.35] 3.27 0.09 Tolerance 95% CI [0.09, 0.10] [0.09, 0.10]

## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
##   For precise Bayes factors, sampling at least 40,000 posterior
##   samples is recommended.
## Warning in summarize_brms(pa_obj_log_pressure, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
summary_pa_obj %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
is_A 101.90*** 6.39 [ 89.74, 115.49] 1.000 [0.94, 1.07] 0.000 >100 Overwhelming Evidence 1.001 7210 14382
is_B 137.60*** 7.65 [123.21, 153.70] 1.000 [0.94, 1.07] 0.000 >100 Overwhelming Evidence 1.000 10736 18368
is_A:pressure_self_cw 0.99 0.05 [ 0.87, 1.09] 0.609 [0.94, 1.07] 0.764 0.006 Very Strong Evidence for Null 1.000 35375 26044
is_B:pressure_self_cw 0.99 0.04 [ 0.90, 1.08] 0.616 [0.94, 1.07] 0.843 0.005 Very Strong Evidence for Null 1.000 54399 31354
is_A:pressure_partner_cw 1.01 0.04 [ 0.94, 1.09] 0.603 [0.94, 1.07] 0.895 0.004 Very Strong Evidence for Null 1.000 56778 29304
is_B:pressure_partner_cw 1.02 0.06 [ 0.92, 1.17] 0.667 [0.94, 1.07] 0.709 0.006 Very Strong Evidence for Null 1.000 37773 26375
is_A:pressure_self_cb 0.84 0.28 [ 0.43, 1.67] 0.697 [0.94, 1.07] 0.135 0.017 Very Strong Evidence for Null 1.000 17538 23698
is_B:pressure_self_cb 0.94 0.22 [ 0.58, 1.51] 0.607 [0.94, 1.07] 0.210 0.012 Very Strong Evidence for Null 1.000 18333 23028
is_A:pressure_partner_cb 1.25 0.33 [ 0.73, 2.14] 0.797 [0.94, 1.07] 0.136 0.022 Very Strong Evidence for Null 1.000 16840 23213
is_B:pressure_partner_cb 1.20 0.36 [ 0.65, 2.21] 0.722 [0.94, 1.07] 0.139 0.014 Very Strong Evidence for Null 1.000 19355 24488
is_A:day 0.95 0.04 [ 0.87, 1.05] 0.842 [0.94, 1.07] 0.642 0.006 Very Strong Evidence for Null 1.000 83189 29167
is_B:day 0.92 0.04 [ 0.84, 1.01] 0.961 [0.94, 1.07] 0.360 0.018 Very Strong Evidence for Null 1.000 90539 28746
is_A:weartime_self_cw 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.07] 1.000 13.519 Strong Evidence 1.000 76844 32076
is_B:weartime_self_cw 1.00 0.00 [ 1.00, 1.00] 0.973 [0.94, 1.07] 1.000 0.025 Very Strong Evidence for Null 1.000 89204 27836
is_A:weartime_self_cb 1.00 0.00 [ 1.00, 1.00] 0.631 [0.94, 1.07] 1.000 0.012 Very Strong Evidence for Null 1.000 9815 16497
is_B:weartime_self_cb 1.00 0.00 [ 1.00, 1.00] 0.725 [0.94, 1.07] 1.000 0.014 Very Strong Evidence for Null 1.000 12846 21776
sd(is_A) 0.34 0.04 [0.27, 0.44] NA NA NA NA NA 1.000 10546 18623
sd(is_B) 0.29 0.04 [0.23, 0.38] NA NA NA NA NA 1.000 12764 19054
sd(is_A:pressure_self_cw) 0.08 0.07 [0.00, 0.25] NA NA NA NA NA 1.000 14433 21148
sd(is_B:pressure_self_cw) 0.05 0.04 [0.00, 0.17] NA NA NA NA NA 1.000 21197 20973
sd(is_A:pressure_partner_cw) 0.04 0.03 [0.00, 0.13] NA NA NA NA NA 1.000 27025 21802
sd(is_B:pressure_partner_cw) 0.07 0.06 [0.00, 0.27] NA NA NA NA NA 1.000 19755 22073
sigma 0.56 0.01 [0.54, 0.57] NA NA NA NA NA 1.000 76294 29458

Pushing

formula <- bf(
  pa_obj ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    barriers_self_cw + barriers_self_cb + 
    is_A:day + is_B:day + 
    
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log_pushing <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_pushing", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
  check_brms(pa_obj_log_pushing, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log_pushing, integer = TRUE, outliers_type = 'bootstrap')
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                     Term  VIF   VIF 95% CI Increased SE Tolerance
##                     is_A 2.00 [1.88, 2.13]         1.41      0.50
##                     is_B 1.02 [1.00, 1.21]         1.01      0.98
##         barriers_self_cw 1.07 [1.03, 1.13]         1.03      0.94
##         barriers_self_cb 1.07 [1.04, 1.14]         1.03      0.93
##     is_A:pushing_self_cw 1.05 [1.02, 1.12]         1.02      0.96
##     is_B:pushing_self_cw 1.05 [1.02, 1.12]         1.02      0.96
##  is_A:pushing_partner_cw 1.10 [1.07, 1.17]         1.05      0.91
##  is_B:pushing_partner_cw 1.96 [1.85, 2.09]         1.40      0.51
##     is_A:pushing_self_cb 1.83 [1.73, 1.94]         1.35      0.55
##     is_B:pushing_self_cb 1.96 [1.85, 2.09]         1.40      0.51
##  is_A:pushing_partner_cb 1.89 [1.78, 2.01]         1.37      0.53
##  Tolerance 95% CI
##      [0.47, 0.53]
##      [0.83, 1.00]
##      [0.88, 0.97]
##      [0.88, 0.96]
##      [0.89, 0.98]
##      [0.89, 0.98]
##      [0.86, 0.94]
##      [0.48, 0.54]
##      [0.51, 0.58]
##      [0.48, 0.54]
##      [0.50, 0.56]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         88%
##     lognormal          9%
##       weibull          3%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##                   lognormal         72%
##                     tweedie         16%
##  neg. binomial (zero-infl.)          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 13, observations = 1594, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000925345 0.005646173
## sample estimates:
## outlier frequency (expected: 0.00331869510664994 ) 
##                                        0.008155583
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_pushing$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log_pushing, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log_pushing, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
summary_pa_obj %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
is_A 106.40*** 7.11 [ 93.02, 121.57] 1.000 [0.94, 1.06] 0.000 >100 Overwhelming Evidence 1.000 15901 24332
is_B 146.91*** 9.89 [128.51, 168.01] 1.000 [0.94, 1.06] 0.000 >100 Overwhelming Evidence 1.000 19659 25323
barriers_self_cw 0.95** 0.02 [ 0.91, 0.99] 0.996 [0.94, 1.06] 0.647 0.177 Moderate Evidence for Null 1.000 91051 28693
barriers_self_cb 0.82* 0.08 [ 0.68, 0.99] 0.980 [0.94, 1.06] 0.070 0.115 Moderate Evidence for Null 1.000 20765 25131
is_A:pushing_self_cw 1.07 0.04 [ 0.99, 1.15] 0.966 [0.94, 1.06] 0.433 0.046 Strong Evidence for Null 1.000 42743 30873
is_B:pushing_self_cw 1.00 0.03 [ 0.94, 1.06] 0.506 [0.94, 1.06] 0.951 0.007 Very Strong Evidence for Null 1.000 42212 30965
is_A:pushing_partner_cw 0.98 0.02 [ 0.94, 1.04] 0.728 [0.94, 1.06] 0.964 0.007 Very Strong Evidence for Null 1.000 67404 32270
is_B:pushing_partner_cw 1.06 0.04 [ 0.99, 1.14] 0.957 [0.94, 1.06] 0.491 0.035 Strong Evidence for Null 1.000 42448 31807
is_A:pushing_self_cb 0.98 0.28 [ 0.55, 1.73] 0.530 [0.94, 1.06] 0.174 0.015 Very Strong Evidence for Null 1.000 16613 23372
is_B:pushing_self_cb 0.98 0.31 [ 0.52, 1.88] 0.525 [0.94, 1.06] 0.154 0.018 Very Strong Evidence for Null 1.000 17222 23611
is_A:pushing_partner_cb 1.16 0.36 [ 0.63, 2.18] 0.688 [0.94, 1.06] 0.140 0.018 Very Strong Evidence for Null 1.000 15659 23690
is_B:pushing_partner_cb 1.42 0.41 [ 0.79, 2.58] 0.884 [0.94, 1.06] 0.080 0.033 Strong Evidence for Null 1.000 17956 24013
is_A:day 1.00 0.07 [ 0.88, 1.14] 0.523 [0.94, 1.06] 0.661 0.005 Very Strong Evidence for Null 1.000 81551 28339
is_B:day 0.89 0.06 [ 0.78, 1.01] 0.969 [0.94, 1.06] 0.186 0.031 Strong Evidence for Null 1.000 80514 29652
is_A:weartime_self_cw 1.00** 0.00 [ 1.00, 1.00] 0.997 [0.94, 1.06] 1.000 0.237 Moderate Evidence for Null 1.000 74748 33237
is_B:weartime_self_cw 1.00* 0.00 [ 1.00, 1.00] 0.993 [0.94, 1.06] 1.000 0.112 Moderate Evidence for Null 1.000 83919 30123
is_A:weartime_self_cb 1.00 0.00 [ 1.00, 1.00] 0.703 [0.94, 1.06] 1.000 0.014 Very Strong Evidence for Null 1.000 16870 23704
is_B:weartime_self_cb 1.00 0.00 [ 1.00, 1.00] 0.882 [0.94, 1.06] 1.000 0.028 Very Strong Evidence for Null 1.000 19627 25505
sd(is_A) 0.33 0.05 [0.25, 0.44] NA NA NA NA NA 1.000 18055 25973
sd(is_B) 0.34 0.05 [0.26, 0.45] NA NA NA NA NA 1.000 17389 24079
sd(is_A:pushing_self_cw) 0.10 0.05 [0.01, 0.21] NA NA NA NA NA 1.000 9624 11024
sd(is_B:pushing_self_cw) 0.09 0.04 [0.01, 0.17] NA NA NA NA NA 1.000 11822 9435
sd(is_A:pushing_partner_cw) 0.03 0.03 [0.00, 0.10] NA NA NA NA NA 1.000 23687 23076
sd(is_B:pushing_partner_cw) 0.10 0.04 [0.01, 0.18] NA NA NA NA NA 1.001 16684 13022
sigma 0.52 0.01 [0.50, 0.54] NA NA NA NA NA 1.000 59441 29518

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Persuasion

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  
  , brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss_persuasion <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_persuasion", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
  check_brms(mood_gauss_persuasion, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss_persuasion, integer = FALSE)
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                        Term  VIF   VIF 95% CI Increased SE Tolerance
##                        is_A 1.17 [1.13, 1.21]         1.08      0.86
##                        is_B 1.04 [1.02, 1.09]         1.02      0.96
##     is_A:persuasion_self_cw 1.06 [1.04, 1.11]         1.03      0.94
##     is_B:persuasion_self_cw 1.07 [1.04, 1.11]         1.03      0.93
##  is_A:persuasion_partner_cw 1.04 [1.02, 1.09]         1.02      0.96
##  is_B:persuasion_partner_cw 2.85 [2.71, 2.99]         1.69      0.35
##     is_A:persuasion_self_cb 2.80 [2.67, 2.94]         1.67      0.36
##  Tolerance 95% CI
##      [0.82, 0.88]
##      [0.92, 0.98]
##      [0.90, 0.96]
##      [0.90, 0.96]
##      [0.92, 0.98]
##      [0.33, 0.37]
##      [0.34, 0.38]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         47%
##        normal         34%
##   exponential          6%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##        tweedie         41%
##  beta-binomial         38%
##    half-cauchy          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
##  DHARMa outlier test based on exact binomial test with
##  approximate expectations
## 
## data:  model.check
## outliers at both margin(s) = 39, observations = 3688, p-value =
## 2.231e-16
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.007530249 0.014428062
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                           0.01057484
summary_mood <- summarize_brms(
  mood_gauss_persuasion, 
  stats_to_report = stats_to_report,
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = F) 
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
summary_mood %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
is_A 3.60*** 0.12 [ 3.35, 3.84] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.002 3494 6838
is_B 3.88*** 0.14 [ 3.60, 4.15] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.001 4288 9553
is_A:persuasion_self_cw 0.01 0.02 [-0.04, 0.05] 0.633 [-0.11, 0.11] 1.000 0.003 Very Strong Evidence for Null 1.000 53463 30917
is_B:persuasion_self_cw -0.01 0.02 [-0.06, 0.03] 0.717 [-0.11, 0.11] 1.000 0.003 Very Strong Evidence for Null 1.000 47802 29543
is_A:persuasion_partner_cw 0.05* 0.02 [ 0.00, 0.09] 0.982 [-0.11, 0.11] 0.998 0.027 Very Strong Evidence for Null 1.000 45958 31085
is_B:persuasion_partner_cw 0.01 0.02 [-0.04, 0.05] 0.626 [-0.11, 0.11] 1.000 0.003 Very Strong Evidence for Null 1.000 55921 29328
is_A:persuasion_self_cb -0.29 0.33 [-0.97, 0.38] 0.810 [-0.11, 0.11] 0.185 0.023 Very Strong Evidence for Null 1.000 5518 10646
is_B:persuasion_self_cb 0.37 0.33 [-0.30, 1.05] 0.868 [-0.11, 0.11] 0.148 0.033 Strong Evidence for Null 1.001 6266 11337
is_A:persuasion_partner_cb 0.41 0.29 [-0.18, 0.99] 0.916 [-0.11, 0.11] 0.120 0.041 Strong Evidence for Null 1.000 5913 10851
is_B:persuasion_partner_cb -0.20 0.39 [-0.98, 0.57] 0.699 [-0.11, 0.11] 0.205 0.020 Very Strong Evidence for Null 1.000 6005 11939
is_A:day 0.45*** 0.07 [ 0.30, 0.58] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 51520 30176
is_B:day 0.08 0.07 [-0.06, 0.22] 0.871 [-0.11, 0.11] 0.673 0.006 Very Strong Evidence for Null 1.000 54047 28585
sd(is_A) 0.70 0.09 [0.55, 0.90] NA NA NA NA NA 1.000 7045 14067
sd(is_B) 0.79 0.10 [0.63, 1.03] NA NA NA NA NA 1.000 7522 12132
sd(is_A:pushing_self_cw) 0.11 0.06 [0.01, 0.23] NA NA NA NA NA 1.000 10529 10208
sd(is_B:pushing_self_cw) 0.06 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 13285 16387
sd(is_A:pushing_partner_cw) 0.13 0.05 [0.04, 0.24] NA NA NA NA NA 1.000 17381 12874
sd(is_B:pushing_partner_cw) 0.12 0.06 [0.01, 0.25] NA NA NA NA NA 1.000 10329 7714
sigma 0.87 0.01 [0.85, 0.89] NA NA NA NA NA 1.000 52078 28947

Pressure

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  
  , brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss_pressure <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_pressure", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
  check_brms(mood_gauss_pressure, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss_pressure, integer = FALSE)
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                      Term  VIF   VIF 95% CI Increased SE Tolerance
##                      is_A 1.19 [1.15, 1.23]         1.09      0.84
##                      is_B 1.01 [1.00, 1.21]         1.01      0.99
##     is_A:pressure_self_cw 1.01 [1.00, 1.12]         1.01      0.99
##     is_B:pressure_self_cw 1.01 [1.00, 1.12]         1.01      0.99
##  is_A:pressure_partner_cw 1.01 [1.00, 1.42]         1.00      0.99
##  Tolerance 95% CI
##      [0.81, 0.87]
##      [0.83, 1.00]
##      [0.89, 1.00]
##      [0.89, 1.00]
##      [0.70, 1.00]
## 
## Moderate Correlation
## 
##                      Term  VIF   VIF 95% CI Increased SE Tolerance
##  is_B:pressure_partner_cw 8.40 [7.94, 8.89]         2.90      0.12
##     is_A:pressure_self_cb 8.05 [7.61, 8.52]         2.84      0.12
##  Tolerance 95% CI
##      [0.11, 0.13]
##      [0.12, 0.13]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         47%
##        normal         34%
##   exponential          6%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##        tweedie         41%
##  beta-binomial         38%
##    half-cauchy          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
##  DHARMa outlier test based on exact binomial test with
##  approximate expectations
## 
## data:  model.check
## outliers at both margin(s) = 41, observations = 3684, p-value <
## 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.007998021 0.015068026
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                           0.01112921
summary_mood <- summarize_brms(
  mood_gauss_pressure, 
  stats_to_report = stats_to_report,
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = F) 
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
summary_mood %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
is_A 3.64*** 0.13 [ 3.37, 3.89] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.001 2966 6814
is_B 3.88*** 0.14 [ 3.60, 4.17] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.003 3288 6776
is_A:pressure_self_cw -0.06 0.06 [-0.17, 0.06] 0.835 [-0.11, 0.11] 0.845 0.005 Very Strong Evidence for Null 1.000 58461 29843
is_B:pressure_self_cw 0.01 0.06 [-0.10, 0.12] 0.542 [-0.11, 0.11] 0.957 0.003 Very Strong Evidence for Null 1.000 44197 29836
is_A:pressure_partner_cw 0.05 0.06 [-0.06, 0.16] 0.808 [-0.11, 0.11] 0.864 0.005 Very Strong Evidence for Null 1.000 42613 30902
is_B:pressure_partner_cw 0.04 0.06 [-0.07, 0.15] 0.754 [-0.11, 0.11] 0.907 0.004 Very Strong Evidence for Null 1.000 68143 31795
is_A:pressure_self_cb -0.16 0.64 [-1.48, 1.11] 0.600 [-0.11, 0.11] 0.142 0.015 Very Strong Evidence for Null 1.000 10697 15995
is_B:pressure_self_cb 0.12 0.57 [-1.03, 1.31] 0.588 [-0.11, 0.11] 0.157 0.017 Very Strong Evidence for Null 1.000 10971 16363
is_A:pressure_partner_cb 0.22 0.50 [-0.76, 1.27] 0.674 [-0.11, 0.11] 0.164 0.016 Very Strong Evidence for Null 1.000 10831 15307
is_B:pressure_partner_cb -0.12 0.72 [-1.61, 1.35] 0.566 [-0.11, 0.11] 0.126 0.016 Very Strong Evidence for Null 1.000 11345 16597
is_A:day 0.41*** 0.07 [ 0.28, 0.55] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 49901 30936
is_B:day 0.09 0.07 [-0.05, 0.22] 0.897 [-0.11, 0.11] 0.647 0.006 Very Strong Evidence for Null 1.000 50514 29576
sd(is_A) 0.71 0.09 [0.57, 0.93] NA NA NA NA NA 1.001 5363 10000
sd(is_B) 0.81 0.10 [0.65, 1.06] NA NA NA NA NA 1.001 6354 13168
sd(is_A:pushing_self_cw) 0.12 0.06 [0.01, 0.25] NA NA NA NA NA 1.001 9570 7098
sd(is_B:pushing_self_cw) 0.07 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 12893 14758
sd(is_A:pushing_partner_cw) 0.13 0.05 [0.04, 0.24] NA NA NA NA NA 1.000 13798 9709
sd(is_B:pushing_partner_cw) 0.12 0.06 [0.02, 0.25] NA NA NA NA NA 1.000 12399 9266
sigma 0.87 0.01 [0.85, 0.89] NA NA NA NA NA 1.000 43311 29442

pushing

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    barriers_self_cw + barriers_self_cb + 
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  
  , brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss_pushing <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_pushing", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
  check_brms(mood_gauss_pushing, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss_pushing, integer = FALSE)
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                     Term  VIF   VIF 95% CI Increased SE Tolerance
##                     is_A 1.35 [1.29, 1.43]         1.16      0.74
##                     is_B 1.01 [1.00, 1.48]         1.01      0.99
##         barriers_self_cw 1.09 [1.05, 1.15]         1.04      0.92
##         barriers_self_cb 1.09 [1.05, 1.15]         1.04      0.92
##     is_A:pushing_self_cw 1.05 [1.02, 1.12]         1.02      0.95
##     is_B:pushing_self_cw 1.08 [1.05, 1.14]         1.04      0.93
##  is_A:pushing_partner_cw 1.09 [1.05, 1.15]         1.04      0.92
##  is_B:pushing_partner_cw 1.50 [1.42, 1.58]         1.22      0.67
##     is_A:pushing_self_cb 1.54 [1.46, 1.62]         1.24      0.65
##  Tolerance 95% CI
##      [0.70, 0.77]
##      [0.68, 1.00]
##      [0.87, 0.95]
##      [0.87, 0.95]
##      [0.89, 0.98]
##      [0.88, 0.96]
##      [0.87, 0.95]
##      [0.63, 0.70]
##      [0.62, 0.68]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         41%
##        normal         41%
##   exponential          6%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##        tweedie         41%
##  beta-binomial         38%
##    half-cauchy          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
##  DHARMa outlier test based on exact binomial test with
##  approximate expectations
## 
## data:  model.check
## outliers at both margin(s) = 11, observations = 1776, p-value =
## 0.00112
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.003095804 0.011055146
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.006193694
summary_mood <- summarize_brms(
  mood_gauss_pushing, 
  stats_to_report = stats_to_report,
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = F) 
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
summary_mood %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
is_A 3.66*** 0.12 [ 3.41, 3.90] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 7565 14933
is_B 3.84*** 0.14 [ 3.55, 4.13] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 7731 13693
barriers_self_cw -0.29*** 0.03 [-0.35, -0.23] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 78719 29655
barriers_self_cb -0.50* 0.22 [-0.95, -0.05] 0.984 [-0.11, 0.11] 0.039 0.167 Moderate Evidence for Null 1.000 8107 14716
is_A:pushing_self_cw 0.09 0.05 [ 0.00, 0.18] 0.969 [-0.11, 0.11] 0.659 0.034 Strong Evidence for Null 1.000 42568 31360
is_B:pushing_self_cw -0.01 0.04 [-0.09, 0.06] 0.640 [-0.11, 0.11] 0.990 0.005 Very Strong Evidence for Null 1.000 48493 32295
is_A:pushing_partner_cw 0.09 0.04 [ 0.00, 0.18] 0.974 [-0.11, 0.11] 0.678 0.039 Strong Evidence for Null 1.000 32188 29785
is_B:pushing_partner_cw 0.13* 0.05 [ 0.02, 0.23] 0.985 [-0.11, 0.11] 0.376 0.077 Strong Evidence for Null 1.000 33217 29086
is_A:pushing_self_cb -0.23 0.48 [-1.20, 0.73] 0.681 [-0.11, 0.11] 0.162 0.015 Very Strong Evidence for Null 1.000 11236 17666
is_B:pushing_self_cb 0.90 0.67 [-0.47, 2.27] 0.905 [-0.11, 0.11] 0.053 0.047 Strong Evidence for Null 1.000 9401 14271
is_A:pushing_partner_cb 0.77 0.52 [-0.29, 1.83] 0.927 [-0.11, 0.11] 0.056 0.042 Strong Evidence for Null 1.001 9141 15831
is_B:pushing_partner_cb -0.32 0.63 [-1.57, 0.96] 0.693 [-0.11, 0.11] 0.124 0.022 Very Strong Evidence for Null 1.000 11283 17409
is_A:day 0.34*** 0.10 [ 0.15, 0.54] 1.000 [-0.11, 0.11] 0.008 1.941 Weak Evidence 1.000 70757 30332
is_B:day 0.08 0.10 [-0.11, 0.27] 0.793 [-0.11, 0.11] 0.609 0.006 Very Strong Evidence for Null 1.000 69006 28830
sd(is_A) 0.64 0.09 [0.50, 0.84] NA NA NA NA NA 1.001 10062 18475
sd(is_B) 0.79 0.10 [0.63, 1.04] NA NA NA NA NA 1.000 10080 17533
sd(is_A:pushing_self_cw) 0.09 0.07 [0.01, 0.24] NA NA NA NA NA 1.000 12455 15086
sd(is_B:pushing_self_cw) 0.07 0.05 [0.00, 0.20] NA NA NA NA NA 1.000 14139 16967
sd(is_A:pushing_partner_cw) 0.12 0.06 [0.01, 0.24] NA NA NA NA NA 1.000 13704 10047
sd(is_B:pushing_partner_cw) 0.14 0.07 [0.02, 0.29] NA NA NA NA NA 1.000 14270 12569
sigma 0.83 0.01 [0.80, 0.86] NA NA NA NA NA 1.000 60197 29169

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}

df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}

Persuasion

formula <- bf(
  is_reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")

  , brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance_persuasion <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_persuasion", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
  check_brms(is_reactance_persuasion)
  DHARMa.check_brms.all(is_reactance_persuasion, integer = FALSE)
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                        Term  VIF   VIF 95% CI Increased SE Tolerance
##                        is_A 1.82 [1.71, 1.96]         1.35      0.55
##                        is_B 1.39 [1.32, 1.49]         1.18      0.72
##     is_A:persuasion_self_cw 1.15 [1.10, 1.23]         1.07      0.87
##     is_B:persuasion_self_cw 1.04 [1.01, 1.14]         1.02      0.96
##  is_A:persuasion_partner_cw 1.06 [1.02, 1.14]         1.03      0.95
##  is_B:persuasion_partner_cw 1.51 [1.43, 1.62]         1.23      0.66
##     is_A:persuasion_self_cb 1.78 [1.67, 1.91]         1.33      0.56
##  Tolerance 95% CI
##      [0.51, 0.58]
##      [0.67, 0.76]
##      [0.82, 0.91]
##      [0.88, 0.99]
##      [0.88, 0.98]
##      [0.62, 0.70]
##      [0.52, 0.60]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        normal         44%
##        cauchy         12%
##         gamma          9%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##      bernoulli         81%
##  beta-binomial         12%
##       binomial          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 5 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with
##  approximate expectations
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0000334886 0.0073476538
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.001322751
summary_is_reactance <- summarize_brms(
  is_reactance_persuasion, 
  stats_to_report = stats_to_report,
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
summary_is_reactance %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
is_A 0.18*** 0.08 [0.07, 0.43] 1.000 [0.83, 1.20] 0.000 36.028 Very Strong Evidence 1.000 24897 29836
is_B 0.30** 0.12 [0.12, 0.65] 0.999 [0.83, 1.20] 0.004 2.802 Weak Evidence 1.000 29961 29498
is_A:persuasion_self_cw 0.90 0.13 [0.64, 1.18] 0.773 [0.83, 1.20] 0.669 0.071 Strong Evidence for Null 1.000 22239 26651
is_B:persuasion_self_cw 0.78 0.12 [0.54, 1.02] 0.965 [0.83, 1.20] 0.310 0.293 Moderate Evidence for Null 1.000 23382 19058
is_A:persuasion_partner_cw 0.89 0.16 [0.58, 1.29] 0.740 [0.83, 1.20] 0.579 0.077 Strong Evidence for Null 1.000 30521 22671
is_B:persuasion_partner_cw 1.27 0.24 [0.86, 1.97] 0.897 [0.83, 1.20] 0.356 0.165 Moderate Evidence for Null 1.000 25245 25126
is_A:persuasion_self_cb 5.13 4.34 [0.96, 35.39] 0.972 [0.83, 1.20] 0.025 0.569 Weak Evidence for Null 1.000 18611 24445
is_B:persuasion_self_cb 2.76 2.14 [0.62, 15.25] 0.910 [0.83, 1.20] 0.078 0.216 Moderate Evidence for Null 1.000 20515 24894
is_A:persuasion_partner_cb 1.58 1.13 [0.40, 7.92] 0.746 [0.83, 1.20] 0.169 0.114 Moderate Evidence for Null 1.000 16564 19332
is_B:persuasion_partner_cb 2.13 1.94 [0.34, 14.06] 0.798 [0.83, 1.20] 0.109 0.122 Moderate Evidence for Null 1.000 18086 24106
is_A:day 2.57 1.36 [0.91, 7.43] 0.962 [0.83, 1.20] 0.058 0.194 Moderate Evidence for Null 1.000 56235 31190
is_B:day 1.04 0.59 [0.33, 3.17] 0.523 [0.83, 1.20] 0.248 0.045 Strong Evidence for Null 1.000 52331 30861
sd(is_A) 1.38 0.42 [0.68, 2.40] NA NA NA NA NA 1.000 14059 19287
sd(is_B) 1.36 0.33 [0.82, 2.18] NA NA NA NA NA 1.000 16662 23957
sd(is_A:persuasion_self_cw) 0.31 0.18 [0.03, 0.73] NA NA NA NA NA 1.001 9510 11805
sd(is_B:persuasion_self_cw) 0.35 0.26 [0.02, 0.95] NA NA NA NA NA 1.000 5934 13700
sd(is_A:persuasion_partner_cw) 0.47 0.27 [0.05, 1.19] NA NA NA NA NA 1.000 12291 11597
sd(is_B:persuasion_partner_cw) 0.65 0.29 [0.19, 1.40] NA NA NA NA NA 1.000 16040 15521

Pressure

formula <- bf(
  is_reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")

  , brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance_pressure <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_pressure", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
  check_brms(is_reactance_pressure)
  DHARMa.check_brms.all(is_reactance_pressure, integer = FALSE)
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                      Term  VIF    VIF 95% CI Increased SE Tolerance
##                      is_A 1.83 [1.72,  1.96]         1.35      0.55
##                      is_B 1.02 [1.00,  1.23]         1.01      0.98
##     is_A:pressure_self_cw 1.01 [1.00,  1.59]         1.01      0.99
##     is_B:pressure_self_cw 1.01 [1.00,  2.01]         1.01      0.99
##  is_A:pressure_partner_cw 1.01 [1.00, 36.14]         1.00      0.99
##  is_B:pressure_partner_cw 2.51 [2.33,  2.71]         1.58      0.40
##     is_A:pressure_self_cb 2.97 [2.76,  3.22]         1.72      0.34
##  Tolerance 95% CI
##      [0.51, 0.58]
##      [0.82, 1.00]
##      [0.63, 1.00]
##      [0.50, 1.00]
##      [0.03, 1.00]
##      [0.37, 0.43]
##      [0.31, 0.36]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        normal         34%
##          beta         16%
##        cauchy         12%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##      bernoulli         81%
##  beta-binomial         12%
##       binomial          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 23 observations with a pareto_k > 0.7 in model 'model'.
## We recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with
##  approximate expectations
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0000334886 0.0073476538
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.001322751
summary_is_reactance <- summarize_brms(
  is_reactance_pressure, 
  stats_to_report = stats_to_report,
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
summary_is_reactance %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
is_A 0.15*** 0.06 [0.06, 0.32] 1.000 [0.83, 1.20] 0.000 >100 Overwhelming Evidence 1.000 27201 29563
is_B 0.21*** 0.07 [0.10, 0.40] 1.000 [0.83, 1.20] 0.000 >100 Overwhelming Evidence 1.000 33935 30484
is_A:pressure_self_cw 2.48 1.27 [0.82, 11.13] 0.957 [0.83, 1.20] 0.046 0.724 Weak Evidence for Null 1.000 21167 16301
is_B:pressure_self_cw 1.84 0.67 [0.78, 4.37] 0.935 [0.83, 1.20] 0.096 0.359 Weak Evidence for Null 1.000 27629 24426
is_A:pressure_partner_cw 0.87 0.65 [0.09, 3.62] 0.581 [0.83, 1.20] 0.205 0.128 Moderate Evidence for Null 1.000 24646 22567
is_B:pressure_partner_cw 1.56 0.77 [0.39, 5.39] 0.805 [0.83, 1.20] 0.170 0.128 Moderate Evidence for Null 1.000 30905 28417
is_A:pressure_self_cb 16.91 28.48 [0.78, 1025.11] 0.965 [0.83, 1.20] 0.018 0.537 Weak Evidence for Null 1.000 26287 20773
is_B:pressure_self_cb 23.56* 34.63 [1.93, 874.04] 0.994 [0.83, 1.20] 0.006 1.817 Weak Evidence 1.000 25130 23908
is_A:pressure_partner_cb 1.05 1.44 [0.05, 17.51] 0.513 [0.83, 1.20] 0.105 0.120 Moderate Evidence for Null 1.000 24236 22816
is_B:pressure_partner_cb 0.36 0.65 [0.00, 9.70] 0.725 [0.83, 1.20] 0.071 0.105 Moderate Evidence for Null 1.000 29760 25453
is_A:day 2.24 1.18 [0.79, 6.37] 0.936 [0.83, 1.20] 0.087 0.138 Moderate Evidence for Null 1.000 59228 29462
is_B:day 1.34 0.70 [0.48, 3.72] 0.710 [0.83, 1.20] 0.234 0.050 Strong Evidence for Null 1.000 64117 29915
sd(is_A) 1.46 0.39 [0.84, 2.44] NA NA NA NA NA 1.000 14625 22279
sd(is_B) 1.14 0.28 [0.67, 1.83] NA NA NA NA NA 1.000 15965 24630
sd(is_A:pressure_self_cw) 1.30 1.07 [0.06, 3.91] NA NA NA NA NA 1.000 9827 16955
sd(is_B:pressure_self_cw) 1.03 0.55 [0.17, 2.56] NA NA NA NA NA 1.000 13715 12493
sd(is_A:pressure_partner_cw) 1.81 1.16 [0.13, 4.53] NA NA NA NA NA 1.000 13418 11212
sd(is_B:pressure_partner_cw) 0.74 0.66 [0.03, 2.89] NA NA NA NA NA 1.000 22546 21566

pushing

formula <- bf(
  is_reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    #is_A:plan_self + is_B:plan_self + 
    #is_A:plan_partner + is_B:plan_partner +
    barriers_self_cw + barriers_self_cb + 
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")

  , brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
  , brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance_pushing <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_pushing", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
  check_brms(is_reactance_pushing)
  DHARMa.check_brms.all(is_reactance_pushing, integer = FALSE)
}
## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                     Term  VIF   VIF 95% CI Increased SE Tolerance
##                     is_A 1.85 [1.75, 1.97]         1.36      0.54
##                     is_B 1.06 [1.03, 1.13]         1.03      0.94
##         barriers_self_cw 1.08 [1.05, 1.15]         1.04      0.92
##         barriers_self_cb 1.24 [1.19, 1.31]         1.11      0.81
##     is_A:pushing_self_cw 1.14 [1.10, 1.21]         1.07      0.87
##     is_B:pushing_self_cw 1.01 [1.00, 1.89]         1.00      0.99
##  is_A:pushing_partner_cw 1.05 [1.02, 1.12]         1.02      0.96
##  is_B:pushing_partner_cw 1.53 [1.45, 1.62]         1.24      0.65
##     is_A:pushing_self_cb 1.33 [1.27, 1.40]         1.15      0.75
##  Tolerance 95% CI
##      [0.51, 0.57]
##      [0.89, 0.97]
##      [0.87, 0.96]
##      [0.77, 0.84]
##      [0.83, 0.91]
##      [0.53, 1.00]
##      [0.89, 0.98]
##      [0.62, 0.69]
##      [0.71, 0.79]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        normal         44%
##        cauchy         12%
##       tweedie          9%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##      bernoulli         75%
##  beta-binomial         16%
##       binomial          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with
##  approximate expectations
## 
## data:  model.check
## outliers at both margin(s) = 0, observations = 486, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.000000000 0.007561553
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                                    0
summary_is_reactance <- summarize_brms(
  is_reactance_pushing, 
  stats_to_report = stats_to_report,
  #model_rows_fixed = model_rows_fixed,
  #model_rows_random = model_rows_random,
  #model_rownames_fixed = model_rownames_fixed,
  #model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
summary_is_reactance %>%
  print_df(
    #rows_to_pack = rows_to_pack
    )
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
is_A 0.13*** 0.07 [0.04, 0.33] 1.000 [0.83, 1.20] 0.000 75.098 Very Strong Evidence 1.000 21167 22208
is_B 0.16*** 0.07 [0.06, 0.35] 1.000 [0.83, 1.20] 0.000 >100 Overwhelming Evidence 1.000 28688 28514
barriers_self_cw 0.80 0.17 [0.53, 1.19] 0.865 [0.83, 1.20] 0.388 0.098 Strong Evidence for Null 1.000 43333 29118
barriers_self_cb 2.74 2.06 [0.61, 12.70] 0.910 [0.83, 1.20] 0.077 0.221 Moderate Evidence for Null 1.000 20027 24196
is_A:pushing_self_cw 1.35 0.25 [0.93, 2.00] 0.944 [0.83, 1.20] 0.258 0.217 Moderate Evidence for Null 1.000 35407 26611
is_B:pushing_self_cw 1.44 0.29 [0.96, 2.23] 0.964 [0.83, 1.20] 0.167 0.440 Weak Evidence for Null 1.000 25829 25680
is_A:pushing_partner_cw 0.77 0.24 [0.32, 1.36] 0.815 [0.83, 1.20] 0.338 0.129 Moderate Evidence for Null 1.000 19974 15679
is_B:pushing_partner_cw 1.15 0.23 [0.75, 1.74] 0.752 [0.83, 1.20] 0.517 0.079 Strong Evidence for Null 1.000 35235 26529
is_A:pushing_self_cb 38.45* 65.62 [2.09, 2544.92] 0.993 [0.83, 1.20] 0.005 2.260 Weak Evidence 1.000 12493 16241
is_B:pushing_self_cb 2.37 3.65 [0.13, 75.77] 0.717 [0.83, 1.20] 0.082 0.123 Moderate Evidence for Null 1.000 17952 21560
is_A:pushing_partner_cb 0.08 0.16 [0.00, 4.46] 0.898 [0.83, 1.20] 0.032 0.268 Moderate Evidence for Null 1.000 16536 21879
is_B:pushing_partner_cb 1.25 1.87 [0.05, 26.76] 0.561 [0.83, 1.20] 0.097 0.100 Strong Evidence for Null 1.000 19752 22818
is_A:day 1.28 0.94 [0.29, 5.36] 0.633 [0.83, 1.20] 0.184 0.059 Strong Evidence for Null 1.000 39476 29465
is_B:day 1.40 0.92 [0.39, 5.18] 0.696 [0.83, 1.20] 0.194 0.064 Strong Evidence for Null 1.000 47667 31034
sd(is_A) 1.52 0.52 [0.61, 2.78] NA NA NA NA NA 1.000 8591 8207
sd(is_B) 1.26 0.38 [0.61, 2.16] NA NA NA NA NA 1.000 14530 19907
sd(is_A:pushing_self_cw) 0.26 0.23 [0.01, 0.91] NA NA NA NA NA 1.000 10496 17794
sd(is_B:pushing_self_cw) 0.56 0.30 [0.07, 1.32] NA NA NA NA NA 1.001 8769 8991
sd(is_A:pushing_partner_cw) 0.64 0.46 [0.04, 1.94] NA NA NA NA NA 1.000 10726 13385
sd(is_B:pushing_partner_cw) 0.28 0.25 [0.01, 1.05] NA NA NA NA NA 1.000 14535 18263
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)

report::report_packages()
  • rmarkdown (version 2.29; Allaire J et al., 2024)
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • posterior (version 1.6.0; Bürkner P et al., 2024)
  • extrafont (version 0.19; Chang W, 2023)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • rlang (version 1.1.4; Henry L, Wickham H, 2024)
  • tidybayes (version 3.0.7; Kay M, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • patchwork (version 1.3.0; Pedersen T, 2024)
  • R (version 4.4.2; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • plotly (version 4.10.4; Sievert C, 2020)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • ggridges (version 0.5.6; Wilke C, 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()